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Abstract 

We  present  particle  dynamics  simulations  for  the  response  of  magnetorheological 
(MR)  fluids  upon  application  of  a  magnetic  field.  The  particles  motion  is  considered  to 
be  governed  by  magnetic,  hydrodynamic  and  repulsive  interactions.  Fluid-particle  in¬ 
teractions  are  accounted  for  via  Stokes’  drag  while  inter-particle  repulsions  are  modeled 
through  approximate  hard  sphere  rejections.  In  accordance  with  their  greater  signifi¬ 
cance,  on  the  other  hand,  (linear)  magnetic  interactions  are  fully  simulated.  The  time 
evolution  is  considered  to  be  magnetically  quasi-static  and  magnetostatic  forces  are  de¬ 
rived  from  the  solution  of  (steady)  Maxwell’s  equations,  recomputed  at  each  instant  in 
time.  For  this  we  use  a  potential  theoretic  formulation  where  the  boundary  integral 
equations  are  solved  with  a  fast  multipole  approach.  We  show  that  the  resulting  nu¬ 
merical  codes  can  be  effectively  used  to  study  a  number  of  experimental  observables 
such  as  effective  magnetic  permeabilities  and  response  time-scales  which  are  of  crucial 
importance  in  the  design  of  MR  fluids. 
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1  Introduction 

Magnetorheological  (MR)  and  electrorheological  (ER)  fluids  constitute  important  classes  of 
“smart”  controllable  materials.  The  essential  characteristic  of  MR.  (resp.,  ER.)  fluids  is  that 
they  may  be  reversibly  transformed  from  a.  liquid  state  to  that  of  a.  Bingham  solid  upon 
application  of  a.  magnetic  (resp.,  electric)  field.  Although  the  discovery  of  these  materials 
dates  back  to  over  half  a  century  ago  (Winslow  [1947],  [1949]  for  ER.  &  Rabinow  [1948] 
for  MR),  their  industrial  realization  had,  until  recently,  proved  elusive  due,  in  part,  to 
the  stringent  stability  and  durability  requirements  that  such  applications  demand.  Indeed, 
it  was  only  a.  few  years  ago  that  commercially  viable  fluids  of  this  type  were  achieved 
(Carlson  et.  a.l.  [1990],  [1994]  &  Leventon  [1993])  and  controllable  devices  based  on  this 
technology  are  now  beginning  to  evolve  (see  www.mrfluid.com).  For  the  most  part,  these 
devices  employ  magnetic  fluids  as  these  present  a.  number  of  advantages  over  their  electric 
counterparts.  These  include  higher  achievable  yield  stresses  (20  to  50  times  those  attainable 
with  ER  suspensions),  lower  voltage  demands,  a.  wider  range  of  possible  liquid  carriers  and 
lower  sensitivity  to  contaminants  and  extremes  in  temperature  (Weiss  et  a.l.  [1992],  [1993a], 
[1993b],  Carlson  at  a.l.  [1994],  Jolly  et.  a.l.  [1998a],  [1998b]). 

Typical  ER  and  MR  fluids  consist  of  micron-sized  electrically  conducting  or  magnetically 
permeable  particles  dispersed  in  a.  carrier  oil.  It  is  this  composition  that,  in  fact,  results  in 
their  field-dependent  rheology:  an  applied  field  polarizes  the  particles  which  therefore  move 
to  reduce  the  stored  electromagnetic  energy  of  the  ensemble.  An  energetically  favorable 
arrangement  consists  of  particle  chains  aligned  in  the  direction  of  the  applied  field  and  this, 
in  turn,  gives  rise  to  a.  strong  resistance  to  applied  strains  (on  the  order  of  lOOkPa.  for  MR 
fluids),  see  e.g.  Carlson  et  a.l.  [1990],  Weiss  et  a.l.  [1992],  Jolly  et  a.l.  [1997],  [1998a], 
[1998b],  In  an  attempt  at  understanding  the  main  factors  affecting  this  so-called  “MR. 
(or  ER)-effect”  and  with  a.  view  to  its  possible  enhancement  for  potential  applications,  a. 
myriad  of  experimental  and  theoretical  studies  have  been  carried  out  since  the  effect  was 
first  recognized  (see,  e.g.,  Ta.o  [1992],  Bullough  [1996]  and  references  therein).  The  great 
number  of  variables  that  influence  this  behavior,  however,  conspired  against  the  derivation 
of  an  accurate  quantitative  assessment  of  their  relative  importance  since  experiments  could 
only  explore  small  fractions  of  parameter  space  and  theoretical  work  had  to  rely  on  over¬ 
simplified  models.  As  it  became  practicable,  on  the  other  hand,  it  also  became  clear  that 
numerical  simulations  on  more  complete  models  could  deliver  substantial  information  for 
the  design  of  new  or  improved  fluids.  In  particular,  over  the  last  decade,  a.  number  of 
numerical  studies  were  performed  on  particle  dynamics  models  for  ER  fluids  that  have 
sought  to  account  for  the  different  types  of  interactions  (electric,  hydrodynamic,  Brownian, 
etc)  that  arise  in  the  presence  of  applied  fields  (Klingenberg  et  a.l.  [1989],  Hass  [1993], 
Hess  and  Weider  [1996],  Mohebi  et  a.l.  [1996],  Parthasarathy  and  Klingenberg  [1996]).  In 
this  manner,  the  qualitative  features  of  a.  number  of  experimental  observables  (such  as  the 
path  of  chain  formation)  were  recreated  and  new  insight  was  gained  into  the  mechanisms 
responsible  for  the  ER  effect. 

In  order  to  make  these  models  numerically  tractable,  a.  number  of  simplifications  had  to 
be  performed.  Some  of  them,  such  as  neglecting  Brownian  forces,  could  be  easily  justified 
on  the  basis  of  simple  scaling  arguments  (see  Sec.  2.2).  Others,  however,  were  only  made  for 
computational  convenience.  In  fact,  such  was  the  case,  for  instance,  with  the  inter-particle 
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electromagnetic  forces  as  well  as  with  their  hydrodynamic  behavior.  Although  these  are 
clearly  the  dominant  effects  within  the  system,  their  accurate  mathematical  representation 
required  what  appeared  to  be  a  “ practically  impossible'’'’  ( Hon  riecaze  and  Brady  [p.  2188, 
1991])  solution  of  Stokes  or  Maxwell’s  equations  in  a  highly  oscillatory  geometry  and  were 
therefore,  in  all  of  these  treatments,  replaced  by  point-dipole  and  Stokes  drag  approxima¬ 
tions  respectively.  We  intend  to  show  here,  on  the  other  hand,  that  some  recent  advances  in 
the  development  of  computational  algorithms  do  indeed  allow  for  an  accurate  and  efficient 
treatment  of  such  oscillatory  problems.  While  it  is  still  certainly  true  that  standard  finite- 
element  or  boundary- element  calculations  on  these  models  (where  each  particle,  or  particle 
boundary,  would  comprise  several  elements  and  where  all  of  which  would  contribute  to  the 
overall  held)  are  well  beyond  today’s  computational  capabilities,  we  shall  demonstrate  that 
such  calculations  become  practicable  through  the  use  of  a  fast  multipole  method  (FMM). 

The  basic  idea  behind  FMM  (Rokhlin  [1985]  Greengard  &  Rokhlin  [1987]  and  Green¬ 
gard  [1988])  consists  of  using  multipole  expansions  to  calculate  far-held  contributions  to  a 
slowly  decaying  (e.g.,  Newtonian)  potential  at  any  given  point.  This  observation  can,  in 
fact,  be  effectively  used  in  a  wide  variety  of  (linear)  potential  calculations  by  appealing  to 
the  iterative  solution  of  the  corresponding  boundary  integral  equations  (see  Sec.  3).  For 
the  present  work  we  have  initially  implemented  an  FMM  to  handle  the  full  (linear)  electro¬ 
magnetic  interactions  at  each  time  step  of  a  particle  dynamics  simulation  (Sec.  3),  while  we 
have  kept  with  previous  models  in  treating  the  hydrodynamics  via  Stokes  drag.  We  note 
however,  that  an  analogous  fast  multipole  approach  can  be  used  to  solve  the  full  Stokes 
equations  of  huid  how  (Greengard  &  Kropinski  [1999]).  It  should  also  be  remarked  that,  in 
accordance  with  their  greater  practical  importance,  our  interest  lies  with  MR  huids  which, 
in  fact,  we  use  for  experimental  validation  (Sec.  4).  While  our  model  in  its  present  state 
would,  in  principle,  equally  apply  to  ER  huids,  we  suspect  that  our  assumptions  on  the 
dominance  of  electromagnetic  and  hydrodynamic  forces  constitutes  a  better  approximation 
in  the  case  of  MR  composites,  as  is  evidenced  by  their  aforementioned  enhanced  stability 
properties.  On  the  other  hand,  a  complete  treatment  of  this  magnetic  case  would  neces¬ 
sitate  the  incorporation  of  nonlinear  effects  such  as  magnetic  saturation.  The  means  for 
such  calculations  will  become  clear,  however,  once  we  establish  here  the  possibility  of  fully 
accounting  for  the  linear  magnetic  interactions,  as  these  calculations  could  be  iterated  to 
explore  the  nonlinear  regime.  In  any  case,  these  considerations  and  those  pertaining  to  other 
possible  extensions  (Stokes  equations,  three-dimensional  geometries,  models  of  rheological 
response,  other  higher  order  effects,  etc)  will  be  left  for  future  work. 


2  Equations  of  Motion 


2.1  Governing  Forces  and  Equations 


In  this  section,  we  derive  the  equations  of  motion  for  circular  particles  in  JR2  in  the  presence 
of  an  external  magnetic  held  Ho-  The  motion  of  the  kth  particle  is  described  by  Newton’s 
second  law  of  motion 


Mk 


d2Xk 

dt2 


(2.1) 
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subject  to  the  initial  position  and  velocity 

x°  =  x/,(0)  and  v°  =  v*(0)',  (2.2) 

where  A4k  is  the  mass,  xk  is  the  center,  and  Fk  is  the  total  force  acting  on  the  kth  particle, 
which  occupies  the  region  Qk.  Contributions  to  arise  from  several  sources,  namely  (see 
e.g.,  Klingenberg  et.  al.  [1989],  Parthasarathy  &  Klingenberg  [1996]  and  Mohebi,  Jamasbi 
&  Liu  [1996]): 

•  Magnetic  Forces 

The  magnetic  force  on  Qk  can  be  calculated  from  the  local  field  H  with  the  aid  of 
the  Maxwell  stress  tensor  erMax  =  fikfio[HlRT  —  ||H|2^],  (£=unit  tensor,  fik,  po  =  the 
permeability  for  the  kth  particle  and  the  carrier  oil,  respectively)  as 

F“ag  =  — /  aMax  •  n  dS,  (2.3) 

'2tt/i0  Jgnk 

where  n  is  the  unit  normal  vector  on  dSlk. 


•  Hydrodynamic  Forces 


For  small  Reynolds  number,  the  hydrodynamics  can  be  decoupled  from  the  magneto¬ 
static  problem  and  can  be  approximated  by  the  Stokes’  drag  (see  Batchelor  [1967]  pp. 
244  for  details).  Indeed,  the  hydrodynamic  force  can  be  approximated  using  Oseen’s 
equations  for  flow  due  to  a  moving  body  at  small  Reynolds  number  as 


F "'hydro 
k 


(2.4) 


For  circular  cylinders  for  instance,  the  Stokes  drag  coefficient  D  is  given  by  D  =  1o^7a 
per  unit  length.  Here  rjc  is  the  carrier  oil  absolute  viscosity  and  Re  is  the  Reynolds 
number.  Formulation  (2.4)  has  been  widely  used  to  approximate  the  hydrodynamic 
force  in  dynamic  simulations  for  both  MR  and  ER  fluids  (see,  e.g.,  Klingenberg  et. 
al.  [1989],  Hass  [1993]  and  Mohebi,  Jamasbi  &  Liu  [1996]). 


•  Repelling  Forces 

In  our  simulations,  we  shall  assume  that  both  the  particles  and  the  container  walls 
are  “ hard To  approximate  this  regime,  we  shall  follow  the  work  of  Klingenberg  et. 
al.  [1989]  and  Mohebi  et.  al.  [1996]  and  propose  that  a  “repelling  force”  acts  on  the 
kth  particle  as  it  approaches  others  or  a  wall  of  the  container.  A  simple  model  for 
such  a  force  is  given,  for  instance,  by 

f\ep  =  -iikHlR^ki  e~P\dk‘\  -  nk  pkH*R  e  ^1^'""!,  (2.5) 

i=i 


where  M  is  the  total  number  of  particles,  H o  is  the  intensity  of  the  applied  magnetic 

_  X/  —  xk 

field,  R  is  the  radius  of  the  particle,  /3  >  0  is  the  repelling  parameter,  r ki  =  -rz - r, 

R  -  xk\ 

dki  =  dist(fH-,  0/).  The  wall  repelling  force  uses  n*,,  an  outward  unit  normal  vector 
at  a  point  p  on  the  boundary  of  the  container  SI  where  p  is  closest  to  xk  on  Oil.  and 
|^wall|  _  dist(0j.5  OSl ).  Other  models,  including  polynomial  forms,  have  been  studied 
by  Klingenberg  et.  al.  [1989]. 
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2.2  Dimensional  Analysis 


Our  simulations  will  proceed  on  the  non-dimensionalized  equations  of  motion  which  we  now 
derive.  Let  us  begin  by  denoting  the  dimensionless  variables  with  an  asterisk  and  performing 
the  following  scalings  in  equations  (2.1)-(2.2) 


x  =  fix* , 
t  =  rt* 

F™ag  =  F|Frs’*, 

prep  _  psprep,* 

*  k  —*k*k 


phydro  _ 

■  k  ~~ 


It  then  follows  that 


where 


I)  fi  dSf. 
t  dt* 


„  d2x * 

Gk  7. *2 


dt 


Gk  = 


t  is  the  scaling  time, 

if  =  VkHoR, 

e-/3\du\  _ 

47t  rjc 


F“"rep,* 

k 


—  1^1=1  rkl 


n  k  e 


wall  | 


T  = 


DR 


t<k  log  //„ 


dx 


=  -#  +  F"ag’*  +  F 


rep,* 


(2.6) 


MkR  _  ppR2Hlnk 

~  327r2[log 


is  a  dimensionless  constant  and  pp  is  the  density  of  the  particle  per  unit  length.  For  the 
MR  fluids  we  seek  to  model,  the  particles  are  basically  composed  of  iron  with  diameters  of 
3  to  5  microns  and  pk  —  103  X  pn .  The  carrier  oil  viscosity  r)c  =  0.01  —  0.1  Pa-s  and  the 
applied  magnetic  flux  //0 H o  =  0.001  —  0.1  Tesla,  so  that  Gk  —  10-4.  Thus  the  term  on  the 
left  of  (2.6)  is  several  orders  of  magnitude  smaller  than  those  on  the  right  hand  side,  and 
we  therefore  set  it  to  zero.  Consequently,  the  equations  of  motion  (2. 1 )  (2.2 )  become 


dS?* 

k  _  rimag,1 

dt*  ~  k 


+  F 


rep,* 


and  x°  =  Xfc(O). 


(2.7) 


It  is  to  be  remarked  that  for  MR  fluids  with  characteristics  as  described  above,  the 
thermal  effects  from  the  continuous  phase  molecules  on  the  particles  are  quite  small  for 
rapid  held  application  (Mohebi  et.  al.  [1996]).  Therefore,  the  magnetostatic  forces  largely 
dominate  the  Brownian  forces.  More  precisely,  the  ratio  of  the  Brownian  force  to  the 
magnetic  force  is  approximately  given  by  A  =  ^2T^2 ,  where  K  =  1.38  X  10-23  Joules/K 
is  the  Boltzmann’s  constant  and  T  =  298K  is  the  operating  temperature.  In  our  case,  this 
ratio  is  of  order  10-8  which  justifies  our  neglecting  the  effects  of  Brownian  motion. 


3  Magnetic  Forces 

Clearly,  the  main  challenge  with  the  model  (2.7)  consists  of  the  calculation  of  the  highly 
oscillatory  magnetic  interactions  {F™3®’*}^^.  Indeed,  an  accurate  estimation  of  such  forces 
demands  the  continuous  knowledge  of  the  local  magnetic  field.  H  as  the  particles  rearrange 
themselves,  so  that  Maxwell’s  equations  must  be  resolved  at  each  instant  in  time.  More 
precisely,  let  us  consider  O  C  JR 2  which  is  filled  with  a  non-magnetic  viscous  fluid  and  M 
permeable  circular  particles  Oi,  if  2,  ■  ■  ■ ,  Dm-  Then,  since  the  electromagnetic  time  scale 


6 


Ly-Reitich-Jolly-Banks-Ito 


is  much  shorter  than  that  of  the  motion  itself,  we  may  safely  assume  that  the  fields  are 
governed  by  the  equations  of  magnetostatics,  namely, 

V  ■  B  =  0,  and  V  X  H  =  J,  (3.1) 

where  J  is  the  free  current,  H  is  the  magnetic  field  and  B  is  the  magnetic  induction.  In  our 
case,  the  material  constituents  are  isotropic  and  hence 


B  —  Bj  -f-  jt^gH  —  fi H, 


(3.2) 


where  B;  is  the  intrinsic  induction  caused  by  magnetization,  po  is  the  permeability  of  free 
space  and 


pk  in  the  kth  particle 

po  in  the  carrier  oil. 


(3.3) 


Notice  that  {pk}kLi  are  not  necessarily  identical  and  are  substantially  larger  than  po, 
{[ik  ~  2000/io)-  In  general,  pk  is  defined  as  a  function  of  H  to  model  magnetic  saturation. 
For  moderate  fields,  however,  ///,•  may  be  accurately  approximated  as  a  constant  and  we 
shall  take  this  as  our  working  assumption.  Moreover,  for  the  composites  we  consider,  no  free 
currents  are  present  in  the  domain.  Thus  the  second  equation  in  (3.1)  becomes  V  X  H  =  0 
which,  in  turn,  implies  that  the  magnetic  held  can  be  written  in  terms  of  a  scalar  potential 

fa 

H  =  -Vfa  (3.4) 

As  a  result,  equations  (3. 1 )  (3. 1)  can  be  simplified  to 


V  •  (/iV#  =  0, 


(3.5) 


with  p  given  by  (3.3)  and  pki  po  constant.  Note  that  (3.3)  implies  that  equation  (3.5)  has 
highly  oscillatory  coefficients.  Also,  of  course,  equation  (3.5)  encodes  the  continuity  of  the 
magnetic  potential  (j)  and  of  the  normal  component  of  B.  That  is,  for  any  k  =  1,2,...,  M, 


lina  fap) 

p  — *  Oilk 
p  €  0* 


p  - 
p 


lim 
-  dttk 
€ 


0) 


4(p) 


lim  j ik 
p  —  OQk 

p  €  0* 


d<j> 
9flk,  p 


(P) 


P 


y  9(j> 

hm  po  — 

—  dttk  onkp 


(P) 


p  e 


(3.6) 

(3.7) 


where  nk,v  is  an  outward  unit  normal  vector  at  p  6  dSlk  and  =  is  the  complement 

of  Q,k- 


3.1  Integral  Equation  Formulation  and  Boundary  Element  Discretization 

Although,  as  we  said,  the  coefficients  of  equation  (3.5)  are  rapidly  changing  in  space,  they 
do  remain  constant  in  each  Slk-  Thus  the  overall  potential  can  be  derived  from  appropriate 
charge  densities  supported  on  the  boundaries  of  the  particles.  These  densities  satisfy  certain 
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integral  equations  which  are,  in  principle,  amenable  to  solution  by  finite  (boundary)  element 
approximation.  As  we  shall  discuss  below  (Sec.  3.2),  the  difficulties  associated  with  the  high 
computational  cost  of  classical  boundary  element  approximation  for  this  kind  of  problem 
can,  in  fact,  be  overcome  through  the  implementation  of  the  Fast  Multipole  Method. 

To  derive  the  integral  equations,  let  us  denote  by  <91  hi  the  boundary  of  the  domain  Q 
and  impose  the  following  Dirichlet  boundary  condition  on  dflo 


01 9Si0  -  /■ 


(3.8) 


A  potential  0  satisfying  (3.5)-(3.8)  can  be  represented  by  single-layer  and  double-layer 
integrals,  see,  e.g.,  Colton  &  Kress  [1983]  and  Greengard  &  Moura  [1994],  in  the  form 


0(P)  =  /  G(p,q)r}j(q)d9{q)+  -  (p,q)£(q)ds(q), 

~[jdQ.3  JdQ o  o,9 


p  e  fi.  (3.9) 


Here,  ndn0,q  is  the  outward  unit  normal  vectors  at  q  G  ODq  and  G(p,q)  =  ^log|p  —  q\  is 
the  fundamental  solution  of  the  Laplace’s  equation  in  JR2 .  The  functions  r/j’s  on 
and  £  on  dflo  are  appropriate  (unknown)  surface  densities.  Note  that  the  potential  0  in 
(3.9)  automatically  satisfies 


A0  =  0  on  ilk  for  =  0,1,...,  M. 


and  the  continuity  condition  (3.6)  at  the  interfaces.  In  this  regard,  we  remark  here  that 
our  representation  of  the  potential  0  differs  from  the  standard  one  described  in  most  of  the 
pertinent  literature.  Indeed,  while  a  standard  representation  would  involve  both  single-  and 
double- layer  potentials  at  each  material  interface,  the  formulation  (3.9)  implicitly  guarantees 
the  continuity  of  the  potential  at  particle-fluid  boundaries.  In  addition,  the  use  of  a  sole 
double- layer  integral  on  the  exterior  boundary  ensures  that  the  equations  for  ry? s  and  £  that 
result  from  ( 3.7 )-( 3.8)  constitute  a  system  of  Fredholm  equation  of  the  second  kind.  In  fact, 
using  1  ho  jump  relations  of  potential  theory  (see,  e.g.,  Jaswon  [1977]  and  Colton  &  Kress 
[1983]),  we  obtain  from  equations  ( 3.7 )-( 3.8)  the  following  system  of  integral  equations, 


M 

Tjk{p)  -  2XkJ2  / 

I  Jen. 


d 


dnk 


-G(p,  q)Vj(q)ds(q) 


,p 


-2A  k 


d 


■L 


dG 


£(P)  +  2X!  /  G(p,q)r].j(q)ds(q)  +  2  f 
Jen j  Jd 


dnk}P  JdQ0  dnmoA 
dG 

9fi0  dndUm 


(p,q)t(q)ds(q) 

(p,q)Z(q)ds(q) 


0,  (3.10) 

2  f(p),  (3.11) 


where  Xk  =  -  and  equations  (3.10)-(3.11)  hold  for  p  G  {dDk}^_1  and  8Do,  respec- 

/'/■•  h  /'o 

tively. 

Our  approach  to  the  solution  of  equations  (3. 10)  (3.1  1 )  (at  any  fixed  instant  in  time) 
relies  on  the  inversion  of  their  discretized  version.  To  this  end,  we  divide  each  boundary 
{IT }  w  o,  including  the  exterior  boundary,  into  Nj  disjoint  elements  and  we  denote 
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by  {plj}1^:1  the  midpoint  of  {7]}^  for  j  =  0, 1,  ■  ■  • ,  M  and  N  =  J2y=o  A  j .  To  derive  the 
approximate  equations,  we  assume  that  the  unknown  densities  and  £  in  ( 3. 10)-(  3.11) 

are  constant  over  each  element,  with  collocation  taken  at  the  midpoints.  That  is, 

/,  U(p,q)g(q)ds(q )  A  g(p))  j  ^  U{p,q)ds{  q) 

T?  T? 


for  p  being  either  r]j  or  £  and  U  being  one  of  the  kernels  in  (3. 10)-(3. 1 1) .  Thus,  denoting 
rjj  =  rjj{plj)  and  =  £(pq),  equations  (3. 10)-(3. 1 1)  can  be  transformed  into  a  matrix 
equation 

(I  +  Mjv)Ujv  =  F  N.  (3.12) 

Here,  I  is  the  N  X  N  identity  matrix, 


M  n 


A  B 
C  D  ’ 


Ujv 


0  1  <  i  <  N 

2/(po)  N  +  l<i<N  +  N0, 


where  N  denotes  J2y=i  Aj, 


$  = 

1 

_ 1 

i  Vk  - 

1 

i— 

1 _ 

for  1  <  k  <  M,  £  = 

'  6  ” 

6 

VM 

1 

•  2^ 
Sr 

_ 1 

.  6v0  _ 

A  $(]>*) 

m.pi) 

c  $U) 

Dfti) 


M  N, 


(. Pk  ~  g)  '  nk,p 
’1)  \Pk  ~  <l\‘ 


w  3=1  1=1  J  ^ 

Wo  , 


^k 
7 r 


92 


To  dHk,pidHq 
ttii  i  log  bo  - 

»_1  7  -I  A 


log  \p\  -  q\ds(q), 


i=i 
M  N 


3=1 1=1 
N0 


hAL 


n3 

(q  KAA‘.q 


ft!  H  bo  -  ?l5 


ds(q). 


(3.13) 

(3.14) 

(3.15) 

(3.16) 


3.2  Numerical  Implementation 

To  solve  for  Ujv  in  (3.12),  we  use  an  iterative  method,  namely  GMRES  (Generalized  Mini¬ 
mal  RESidue).  As  an  iterative  solver,  GMRES  demands  the  repeated  calculation  of  products 
(1+  Mjv)Ujv-  The  matrix  Mjv  in  our  linear  system  (3.12)  is  fully  dense  so  that  multiplica¬ 
tion  of  Mjv  and  Ujv  would  require  0(N2)  operations.  However,  as  we  describe  below,  by 
exploiting  the  physical  nature  of  the  underlying  magnetostatic  problem  it  is  possible  to  ac¬ 
curately  approximate  these  matrix  vector  multiplication  by  a  procedure,  the  Fast  Multipole 
Method  (FMM),  whose  operation  count  is  only  of  (9(  A  log  A). 
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The  FMM  algorithm  was  introduced  in  Rokhlin  [1985]  for  rapid  solution  of  integral 
equations  in  potential  theory  and  later  extended  by  Greengard  &  Rokhlin  [1987]  to  fast 
evaluation  of  Coulombic  interactions  in  large-scale  systems  of  particles.  A  typical  calculation 
of  the  field  at  a  given  point  due  to  a  distribution  of  N  charges  (MjyUjv,  in  our  notation) 
can  be  broken  into  two  parts,  the  far  held  and  the  near  held  calculations.  However,  since 
the  Coulombic  potential  decays  rather  slowly  (logarithmically),  the  far-held  calculation 
must,  in  principle,  account  for  all  interactions  among  the  particles  and,  consequently,  the 
number  of  operations  for  calculating  the  held  is  0{N2).  The  goal  of  the  fast  multipole 
method  is  to  reduce  this  number  of  operations  while  accounting  for  all  interactions  and 
maintaining  any  desirable  order  of  accuracy.  To  achieve  this,  the  basic  strategy  consists 
of  clustering  the  charges  at  different  spatial  lengths  to  allow  for  the  computation  of  the 
interactions  between  clusters  by  using  multipole  expansions.  Near-held  interactions  are 
computed  directly.  This  systematic  process  is  the  key  ingredient  in  reducing  the  number  of 
operations  in  the  numerical  matrix- vector  multi  plica  t  ions.  Complete  details  of  the  FMM 
algorithm  can  be  found  in  Greengard  [1988].  For  other  applications  and  implementations 
of  the  FMM  see  also  Rokhlin  [1985]  [1990],  Greengard  &  Rokhlin  [1987]  [1989],  Greengard 
&  Gropp  [1990],  Greengard  &  Strain  [1990],  Engheta  et.  al.  [1992],  Murphy  et.  al.  [1993], 
Greengard  &  Moura  [1994],  Jones  et.  al.  [1994],  McKenney  et.  al.  [1995],  and  Greengard 
&  Kropinski  [1999]. 

3.3  Error  Analysis 

In  this  subsection,  we  present  a  simple  error  analysis  for  the  solution  of  the  boundary 
element  method  discussed  in  the  previous  subsection.  Since  we  employ  the  midpoint  rule 
for  all  the  quadratures,  the  order  of  convergence  is  expected  and  will  be  shown  to  be  at 
least  quadratic.  Let  us  start  with  a  system  of  integral  equations  (I  +  M)U  =  F  where 
U,  F  €  L2(T)  and  M  :  i2(r)  — ■>  i2(r)  is  a  compact  operator  (the  L2{V  [-topology  can,  of 
course,  be  substituted  by  others  depending  on  the  regularity  of  U).  We  assume  M  takes 
the  following  form 

MU(jj)  =  K(p,  q)U(  q)dq. 

Also,  the  discretized  boundary  integral  equations  can  be  written  as  an  A  X  A  system 

(I  +  Mn)Un  =  Fn 

with  the  midpoint  collocated  right-hand  side  [Fjsf]i  =  F(pi )  for  1  <  i  <  N ,  where  the  pfs 
are  the  midpoints  of  the  segments  {r8}^_1,  so  that  [{I  +  M)U](pi)  —  (I  +  =  0.  By 

defining  (PjvG)j  =  U(pi)x{Pt},  we  have 

(I  +  Mn){Un  —  PnU)  =  l'\  Ml  —  MnPnU.  (3.17) 

For  sufficiently  smooth  boundaries  and  M  compact,  it  can  be  shown  that  (see  Kress  [1989]) 
||(J  +  Afjv)-1 \\{p-p}  <  A,  independently  of  N.  (3.18) 

It  follows  from  (3.17)  and  (3.18)  that 


\\Un  ~  PnU\\p  <  C'|AJ|  -  M.\ I }.\t  \  p, 


(3.19) 
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for  N  >  0.  Let  E  =  PnMU  —  M ,v  Pn  U  denote  the  locally  truncated  error  at  p.  Since  we 
are  using  the  midpoint  quadrature  rule,  for  p  away  from  the  midpoints 

N 

K(p,  q ) U ( q ) dq  -  J2K(p,pi)U(pi)\Ti\  \p  <  0(N~2), 
i= 1 


\E(p)\P  =  \Jr 


which  implies 


For  p  —  Pi,  we  compute 


II uN  -  PNU\\p  <  c2k.n~2. 

/  K(pi,  q)U(q)dq  exactly  (see  Greengard  and  Moura  [1994]). 

J\ , 


4  Numerical  Results 

In  this  section,  we  first  describe  the  physical  regime  and  parameters  for  the  dynamic  simu¬ 
lations  of  five  MR  fluids.  We  then  discuss  our  results  and  quantify  the  magnetic  responses 
based  on  their  time  scales  and  effective  permeabilities. 


4.1  Microstructure  Evolutions  in  MR  Fluids 

The  dimensions  for  the  rectangular  domain  that  we  consider  for  dynamic  simulation  are 
Lx  X  Ly  =  0.237mm  X  0.1mm.  We  assume  that  all  particles  are  circular,  initially  centered 
randomly  at  .  We  consider  five  samples  of  MR  fluids,  corresponding  to  5%,  10%, 

15%,  20%  and  30%  volume  fractions.  There  are  170,  341,  511,  682,  and  1024  particles 
for  the  respective  samples.  In  addition,  we  assume  the  following  physical  parameters  in  our 
simulations: 


Radius 

Fluid  Density 
Particle  Density 
Fluid  Permeability 
Particle  Permeability 
Viscosity 
Applied  Field 


R  =  1.5  micrometers 
Po  =  103  kg/rn3 
p  =  7  X  po 
p0  =  47t10 N/A2 
Pk  =  2000  X  po. 
rjc  =  0.1  Pa  s 

Ho  =  7.96 F'4  A/m  (poHo  =  0.10  Tesla). 


The  repelling  parameter  in  (2.5)  and  the  Reynolds  number  are  assumed  to  be  40  and  0.01, 
respectively.  For  the  boundary  conditions  on  the  rectangle,  we  assume  <f>  satisfies 


<?(*:  V) 


HgLy 

0  <  x  <  Lx,  y  -  0 

"oMi  j-) 

ljy 

0  ^  y  ^  Ly ,  x  —  Lx 

0 

0  ^  ^  ^  Lx 5  y  —  Ly 

"oMi  f) 

1jv 

0  <  y  <  Ly ,  x  =  0. 

y 


(4.1) 


The  results  for  the  dynamic  simulations  for  10%,  20%  and  30%  volume  fraction  MR,  fluids 
are  displayed  in  Figures  1-3  respectively. 
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Time  =  51 .29  milliseconds 


Time  =  37.69  milliseconds 


Figure  1:  dynamic  simulations  for  MR.  fluid  with  10%  volume  fraction. 


Based  on  these  results,  we  observe  that,  at  first,  particles  form  short  fragmented  chains 
in  the  direction  of  the  applied  field.  Subsequently  these  short  chains  merge  together  and 
form  longer  chains.  As  time  progresses  further,  these  chain-like  clusters  continue  to  lengthen, 
align  and  approach  a  steady  state.  We  also  observe  that  the  time  for  the  columnar  structures 
to  form  and  the  particle  volume  fraction  (<^)  of  the  sample  are  inversely  related. 


4.2  Quantifications  of  Microstructure  Evolutions 

Recent  experimental  work  has  been  conducted  to  indirectly  measure  the  microstructure 
response  in  MR  fluids  using  real-time  permeability  measurements  (Jolly  et.  al  [1997])  and 
analogously  in  ER  fluids  using  permitivity  measurements  (Blackwood  et.  al.  [1994]).  In 
both  cases,  experimental  data  were  fitted  with  exponential  functions  in  an  attempt  to 
identify  the  time  constants  f  for  microstructure  formation.  A  theoretical  estimate  of  such 
constants  can  be  easily  derived  by  consideration  of  the  time  response  of  a  pair  of  point 
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Time  =  25.51  milliseconds  Time  =  30.31  milliseconds 


Figure  2:  dynamic  simulations  for  MR.  fluid  with  20%  volume  fraction. 


dipoles  within  a  viscous  fluid;  this  results  in  (Jolly  et.  al  [1997]) 

(4- 

with  n  =  5/3.  In  this  subsection,  we  undertake  a  numerical  study  of  these  time  scales. 
Following  the  experimental  procedure,  we  first  derive  values  for  such  constants  from  ef¬ 
fective  permeability  calculations.  Finally,  we  estimate  similar  constants  from  a  different 
macroscopic  measurements,  namely  that  of  the  “ average  kinetic  energy 

We  begin  by  examining  the  evolution  of  the  effective  permeability  of  the  MR.  fluid  as  a 
function  of  time.  As  we  said,  we  do  so  with  the  expectation  that  the  effective  permeability 
reflects  the  microstructure  state  within  the  MR  fluid.  The  definition  and  the  formulation 
for  the  effective  permeability,  which  is  based  on  theory  of  homogenization  (Jikov  [1994]), 
are  derived  in  Ly  et.al.  [1998].  The  effective  permeability  is  defined  as  a  matrix 


U  * 


ClTfc  


y~n  +  C2 


/'ell  = 


1^11  Rl2 
R21  R22 
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Time  =  0  milliseconds 


Time  =  1 .138  milliseconds 


Time  =  1.691  milliseconds 


Time  =  10.87  milliseconds 


Time  =  13.65  milliseconds 


Time  =  15.52  milliseconds 


Figure  3:  dynamic  simulations  for  MR.  fluid  with  30%  volume  fraction. 


where 

Reff  <  V<£(x)  >  =  <  /i(x)  V<^(x)  >,  (4.3) 

and  <  ■  >  denotes  the  spatial  average.  Using  Green’s  identity  and  the  boundary  conditions 
(4.1),  equation  (4.3)  yields  explicitly 


Rl2 


1 


M 


P-oLxL^ 


- 


k=  1 


3fi,. 


<jmxdS 


and 


R22  —  Ro  + 


1 


M 


Ho-fxL 


xo^  - 


k= 1 


3fi,. 


(j)nydS 


Because  the  applied  magnetic  held  is  in  the  vertical  direction,  [122  is  a.  more  relevant  quantity 
and  we  denote  by  fi  :=  Ji 22  the  effective  permeability  reflecting  the  overall  magnetic  response 
of  the  MR  fluid.  In  Figure  4,  we  display  the  effective  permeability  as  a  function  of  time.  We 
remark  that  the  effective  permeabilities  obtained  from  the  dynamic  simulations  are  roughly 


14 


Ly-Reitich- Jolly- B  anks-Ito 


40-50  %  smaller  than  those  obtained  from  numerical  and  physical  experiments  in  Simon  et. 
al  [1998].  One  reason  we  suggest  pertains  to  the  percolation  limit  of  particle  separation.  As 
pointed  out  in  the  study  of  Simon  et.  al  [1998],  the  inter-particle  distance  plays  a  major  role 
in  the  effective  permeability  in  MR  fluids:  the  permeability  increases  as  the  inter- particle 
distance  decreases.  For  the  two-dimensional  case,  it  is  necessary  to  allow  the  particle  gap 
to  approach  1%  of  the  particle  radius  in  order  to  achieve  values  comparable  to  those  of 
experimental  results.  In  our  simulations,  the  gaps  between  particles  were  constrained  to  at 
least  4%  of  the  radius.  One  could  allow  the  particles  in  our  study  to  get  closer  to  1%  of  the 
radius  by  refining  the  boundary  integral  element  which,  of  course,  would  lead  to  prolonged 
computing  time. 


Figure  4:  Least  square  fitting  (line)  on  the  computed  effective  permeability  (star)  and  the 
time  scale  0  (solid  square) 

Figure  4  shows  that  the  values  of  the  effective  permeabilities  elevate  faster  for  samples 
with  higher  concentration  of  iron.  To  identify  the  time-scales  that  correspond  to  this  behav- 
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Magnetic  Time  Response:  t^A  *  VolumeFraction  n+C 


Figure  5:  Time  scale  1 1  obtained  from  effective  permeability  and  its  least-square  fitting  for 
the  power  law  (4.2) 


ior,  we  perform  an  exponential  fit  to  a  function  of  the  form  p(t)  —  A  (1  —  exp{--^-))  +  C  for 
each  of  the  considered  samples.  The  values  of  the  time-scales  f  for  the  respective  volume 
fractions  are  displayed  in  Figure  4  by  the  large  solid  squares.  To  determine  the  dependence 
of  the  time  scale  on  the  volume  fraction  tp,  we  use  equation  (4.2)  to  fit  the  data  in  Figure 
4  to  find  that  n  ~  1  (see  Figure  5).  In  this  regard,  our  experiments  indicate  that  the  value 
of  the  power  n  is  rather  insensitive  to  the  initial  configuration  of  the  system. 

Finally,  motivated  by  the  observation  of  initial  rapid  particle  motion  and  subsequent 
slow  rearrangement,  we  have  monitored  the  evolution  of  the  average  particle  kinetic  energy 
E(t)  as  another  macroscopic  measure  of  microstructure  changes.  For  this  we  define  E{t)  = 
^yVfV(f)2  where  A4  is  the  particle  mass  and  V(f)  is  the  average  velocity, 

vm  =  —  V  -**(*-  Mi 

u  Mh  A/ 

Figure  6  displays  the  average  particle  kinetic  energy  E(f)  and  confirms  that  E(t),  associated 
with  the  motion  of  the  particles,  decays  at  a  rate  that  depends  on  the  volume  fraction.  An 
exponential  fit  in  time  E(t )  =  A  exp(—t/t i)  +  C  as  in  the  previous  case  reveals  that 
{/i }  is  approximately  proportional  to  (p~3-3/3.  This  is  in  remarkable  agreement  with  the 
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x  10 


n  Average  Kinetic  Energy  Approximation:  E(t)=A*exp(-t/t1  )+C 


1.5- 


30%  VF;  A=1.5e-ll;  ^=0.0021;  C=1.6e-13 


§  '  X 

°-5V . 
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x  10  sec 
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a5"  V 

n _ 


0  _,0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09 
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15%  VF;  A=8.7e-12;  ^=0.0061;  C=8.6e-13 
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0  ,  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09 

x  1 0  sec 
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10%  VF;  A=5.3e-12;  ^=0.017;  C=3.1e-13 


0  ,  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09 

x  10  sec 


Figure  6:  Numerical  simulation  kinetic  energy  E(t )  (star),  its  least-square  approximation 
(line)  and  the  time  scale  tp  (solid  square). 


experimental  measurements  Blackwood  et.  al.  [1994]  and  Jolly  et.  al.  [1997]  who  have 
found  n  in  (4.2)  to  be  between  2/3  and  4/3. 

5  Summary  and  Conclusions 

We  have  presented  a  computational  technique  to  perform  particle  dynamics  simulations 
of  MR  fluids  upon  application  of  an  external  magnetic  field  that,  for  the  first  time,  fully 
account  for  all  linear  (long-range)  magnetic  interactions.  To  calculate  these  magnetic  forces 
we  solve  a  (highly  oscillatory)  magnetostatic  problem  at  each  instant  in  the  evolution  by 
appealing  to  a  Fast  Multipole  Method  on  a  boundary  integral  formulation.  Additional 
hydrodynamic  and  repulsive  forces  are  accounted  for  by  Stokes  drag  and  approximate  hard- 
sphere/hard-wall  rejections,  respectively.  We  have  effectively  used  the  resulting  numerical 
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Magnetic  Time  Response:  t^A  *  VolumeFraction  "+C 


Figure  7:  Time  scale  1 1  obtained  from  kinetic  energy  and  its  least-square  fitting  for  the 
power  law  (4.2) 


code  to  study  a  number  of  crucial  experimental  observables  (effective  permeability,  response 
time  scale)  and  have  found  our  results  in  good  agreement  with  experimental  data. 
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